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ABSTRACT 

We study the structural evolution of turbulent molecular clouds under the influence of ionizing 
radiation emitted from a nearby massive star by performing a high resolution parameter study with 
the iVINE code. The temperature is taken to be 10 K or 100 K, the mean number density is either 
100 cm~ 3 or 300 cm~ 3 . Furthermore, the turbulence is varied between Mach 1.5 and Mach 12.5, the 
main driving scale of the turbulence is varied between 1 pc and 8 pc. We vary the ionizing flux by 
an order of magnitude, corresponding to allowing between 0.5% and 5% of the mass in the domain 
to be ionized immediately. In our simulations the ionizing radiation enhances the initial turbulent 
density distribution and thus leads to the formation of pillar-like structures observed adjacent to HII 
regions in a natural way. Gravitational collapse occurs regularly at the tips of the structures. We 
find a clear correlation between the initial state of the turbulent cold cloud and the final morphology 
and physical properties of the structures formed. The most favorable regime for the formation of 
pillars is Mach 4 — 10. Structures and therefore stars only form if the initial density contrast between 
the high density unionized gas and the gas that is going to be ionized is lower than the temperature 
contrast between the hot and the cold gas. The density of the resulting pillars is determined by a 
pressure equilibrium between the hot and the cold gas. A thorough analysis of the simulations shows 
that the complex kinematical and geometrical structure of the formed elongated filaments reflects 
that of observed pillars to an impressive level of detail. In addition, we find that the observed line-of 
sight velocities allow for a distinct determination of different formation mechanisms. Comparing the 
current simulations to previous results and recent observations we conclude that e.g. the pillars of 
creation in M16 formed by the mechanism proposed here and not by the radiation driven implosion 
of pre-existing clumps. 

Subject headings: stars: formation, ISM: structure, turbulence, ultraviolet: ISM, methods: numerical, 
HII regions, ISM: bubbles, ISM: kinematics and dynamics 



1. INTRODUCTION 



Stars are known to form in turbulent, cold molec- 
ular clouds. When massive stars ignite, their UV- 
radiation ionizes and heats the surrounding gas, leading 
to an expanding HII bubble. As soon as the HII re- 
gion breaks through the surface of the molecular cloud 
a low-density, optically thin hole is formed, which re- 
veals the otherwise obscured interior. At the inter- 
face between the HII region and the molecular gas pe- 
culiar structures, often called pillars, are found. The 
most famous examples a re the 'pillars of creation 1 in 
the Eagle Nebula (M16, iHester et ail 1 1 9961 ). There is 
also wide-spread e vidence for star formation at the tips 
of the pillars (e.g. | Sugitani et al.ll20"02t iThompson et al.l 
120021 ISugitani et al.l l2007t lOliveiral 120081) . Since the 
launch of the Spitzer Space Telescope a wealth of 
highly resolved observations of the peculiar, pillar or 
trunk like structures observed around the hot, ion- 
ized HH-regions around massive stars and the star for- 
mation in this tr unks has become available, e.g. in 
the Orion clouds (IStan ke et al.l 120021 iLee fc Chenl 120071: 
Bowler et al.l 120091 ) . the Carina n ebula (ISmith et al.l 
2000D , the Elephant Trunk Nebula (|Reach et al.ll2004D . 
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the Trifi d Nebula (iLefloch et all 12002}) . the Rosette 
Nebu la (ISchneider et al.l l2010l )~]VI16 (lAndersen et al 



2004[) . M17 dJiang et aLll2002D, 30 Dor JWalborn et al 



20021 ) and the SMC (|Gouliermis et~aH l2007f ) 
dition, several r ecent observations of bright 
clouds (BRCs 



In ad- 
rimmed 

(lUrauhart et al.l 120091: IChauhan et"aH 



120091 : iMorgan et all I2009L l2010h have been carried 
out. An interesting aspect is the surprisingly spher- 
ical shape of many observed nebulae, especially in 
RCW 120, 'the perf ect bubble' (peharveng et al.l 120091: 
Zavagn o'et~aIll201J|). Other regi ons, like e.g. RCW 



79 (IZavagno et al.ll2006l). RCW82 (iPomares et al.l l2009h 



RCW 108 (iComeron fc Schneiderl 12007ft a nd Sh 104 
(jDeharveng et al.l 120031 : iRodon et al.1 1201CD share this 
morphology. 

In general, the pillars point like fingers towards the 
ionizing source and show a common head-to-tail struc- 
ture. Most of the mass is concentrated in th e head which 
has a bright rim facing the young stars (e.g. iGahm et al.1 
2006). Thin, elongated pillars connect the head with the 
main body of the molecular cloud. They have typical 
widths of 0.1 — 0.7pc and are 1 — 4 pc long (jGahm et al.1 
120061: ISchuller et al.l 12001 ) . The observations show that 
the pillars are not smooth, b ut show sma ll scale struc- 
ture, filaments and clumps ([Pound! [l998| ). Some fila- 
ments run diagonal across the pillars, suggesting a com- 
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plex twist into a helical structure ([Carlgvist et al.| [2003). 
This is also supported by spectroscopic measurements of 
the line-of-sight (LOS) gas velocity: the pillars show a 
bulk motion away from the ionizing stellar sources with 
a superimposed complex shea r flow that could b e inter- 
preted as corkscrew rotation (jGahm et al.ll2006l ). Occa- 
sionally, close to the tip of the head small spherical gas 
clumps are observed to break off and float into the hot 
HII region. These so called evaporating gaseous globules 
(EGGs) have been found with HST e.g . in the Eagle 
Nebula (jMcCaugh rean fe Andersenl [20021. If stars with 
surrounding gas discs happen to form in these clumps 
they transform into evaporating proto-plane tary discs, so 
called proplyds. More recent observations (|Gahm et al.l 
l20Q7f) have revealed a wealth of even smaller sized glob- 
ules or globulettes, which are in general not bright 
rimmed and show no detectable sign of star formation. 
Direct signatures of star formation are found in the head 
of the pillars, e.g. through jets from obscured proto-stars 
piercing throug h the surface into the HII region (e.g. 
in Eta Carina, ISmith et a.1.1 l2000h . Whether these jets 
are preferentially aligned perpendicular to the trunk is a 
matt er of current debate (jRaga et al.ll2010t ISmith et al.1 
f2010h . 

On the theoretical side, early models of pillar formation 
suggested that they form by Rayleigh- Taylor instabilities 
when the expanding hot, low-d ensity HII regi on radially 
accelerates the cold, dense gas (|Friemanl ll954). This has 
been ruled out b y the observat ions of the complex flows 
inside the pillars (|Poundlll998l ). 

Another possibility is the collect and collapse model. 
Here, the radiation sweeps up a larg e shell, which then 
fragments to form stars and pillars (Elmcgre en fc Ladal 
[l977t Klein et al.lll98Ct ISandford et al.lll982lh However, 
the tim escales (> 5Myr) and masses (> 1000 Mp) in- 
volved (|Elmegreen et al.l 119951 : IWunsch fc Paloul l2001h 
are much larger than in M16. Therefore, this is a more 
likely scenario for supernova-driven shells. 

A third scenario is th e radiation-driven implosion 
(RDI, e.g. iBertoldil 119891) of pre-existing dense cores. 
This has been studied in great detail with numerical 
simul ations fe.g. lLefloch fe Lazareflll 994: Will iams et al.l 
120011 ). More recently iKessel-Devnet fe Burkertl (|2003f ) 
presented three-dimensional RDI simulations with a 
smoothed-particles-hydrodynamics (SPH) code and were 
able to show that an otherwise gravitationally marginally 
stable sphe re can be driven into collap se by ionizing ra- 
diation. In iGritschneder et al.l (j2009aL hereafter G09a) 
we showed that marginally stable density enhancements 
get triggered into for ming stars in cases with high as well 
as low ionizing flux. iMiao et"al"1 (|2009l ) further analyzed 
this RDI-scenario with a SPH-based radiative transfer 
scheme. They show that there is an evolutionary se- 
quence, depe nding on the initia l size o f t he cloud, as 
sugges ted by iLefloch fc Lazarefll (|1994D . iBisbas et al.l 
(2009) studied the implosion of a single clump with a 
new ray-tracing scheme, based on the HEALPix algo- 
rithm. An new implementation in the adaptive-mesh- 
refincmcnt (AMR) code FLASH usi ng the hybrid char- 
acteristics ray tracing was achieved bylPeters et al.l ((2009L 
120101) . Very recently IMackev fc Liml (|2010l ) were able to 



ing recipes in a grid code but are missing the effects of 
self-gravity. 

Recently, the foc us has moved towards the ionization of 
the turbulent ISM. lMellema et al.l (|2006f) reproduced the 
observed morphologies of HII regions by ionizing a tur- 
bulent medium using a gr id code without the inclusion of 
gravity. iDale et all (|2007l) used an SPH code to compare 
the gravitational collapse of a molecular cloud with and 
without ionization. They found slightly enhanced star 
formation in the simulation with ionization. The inclu- 
sion of ionization in a grid cod e in combination with a 
magnetic field was discussed bv lKrumholz et afl {2007). 
A homogenous magnetic field leads to a non-spheric HII- 
region, as the gas is held back by the magnetic field lines 
and an oval shaped bubble develops. However, all these 
simulations focussed on the evolution of the entire HII- 
region and therefore lack the resolution for a detailed 
kinematical and structural analysis of the pillars. For 
a comprehensive quantitative study we investigated the 
radiative ionization of a turbulent molecular cloud from 
a nearby star cluster with so far unprecedented resolu- 
tion by zooming into a subregion of the cloud. First, we 
focussed on the evolution of the power spectrum in the in- 
teraction zone of a turbulen t cloud affected by ionization 
(IGritschneder et al.l [2009H hereafter G09b). Recently, 
Lora et al.1 (j2009) further investigated the ionization of a 
turbulent cloud by the combination of a two-temperature 
equation of state with gravitational forces and transfer 
of ionizing radiation. They produce pillar-like structures 
including cores. The angular momentum of these cores 
is preferentially aligned perpendicular to the direction of 
the ionizing field. 

A di fferent approach was presented bvlNavakshi n"et al.l 
(2009). They focussed on the momentum transfer in 
regimes of high radiation pressure by combining SPH 
with a Monte-Carlo approach. Another main motiva- 
tion of developing software able to treat ionizing radia- 
tion have been investigations of the re-ionization of the 
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early universe (see e.g. Illiev et al.ll2006l:lPawlik fc Schavd 
2008; lAltav et al.ll2 008t Illiev et al.ll2009f and references 
therein) . 

In this work we investigate for the first time the effect 
of different initial conditions on the ionization of turbu- 
lent molecular clouds. We vary the initial temperature, 
the mean number density, the level of turbulence as well 
as the turbulent scale and the ionizing flux. The struc- 
ture of this paper is as follows. In Sj2]we briefly review the 
concept of ionizing radiation, followed by a short sum- 
mary of the iVINE-code. After that we present the set 
of initial conditions for the parameter study. In fj3] the 
outcome of the different simulations is discussed in de- 
tail. |2] is dedicated to the stability and shape of the 
structures in dependance on the initial conditions. A 
close comparison to the observed masses, morphologies 
and line-of-sight (LOS) velocities is done in Sj5] We draw 
the conclusions in $6] 

2. BASIC APPROACH AND INITIAL CONDITIONS 

2.1. Ionizing Radiation 

As soon as an O star is born it ionizes its surroundings 
with its UV-radiation. This leads to an ionized, hot HII- 
region (Ti on ~ 10 4 K). In the beginning the 'R-type' 
ionization front travels with a speed i>r which is larger 
than the sound speed of the hot gas a; on . This phase 



reproduce the density structure of the main pillar in M16 
by setting up several pre-existing dense cores in a trian- 
gular way. They investigate the effects of different cool- 
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ends as soon as ionization is balanced by recombination 
in the HII- region. The ionized volu me Vs, the so called 
Stromgren sphere (Stromgren 1939), is then given by 



J-Ly 



(1) 



Here, Jl y is the flux of ionizing photons of the source, 
which is assumed to be constant and monochromatic, ocb 
is the recombination coefficient, po is the density of the 
pre-existing, homogeneous gas and m p the proton mass. 
At a larger distance from the star the radiation can be 
assumed to impinge onto a surface in a plane-parallel way 
and thus Eq. [1] simplifies to 



«b(po/"i p ) 



(2) 



where xs is the thickness of the ionized region and -FLy 
is the in-falling ionizing flux per unit time and unit area. 

After a hydrodynamical timescale the hot gas reacts to 
its increased temperature and pressure. The pressure of 
an ideal one-atomic gas is given by 



P = p- 



> 

pm p 



(3) 



where p is the density, fee is the Boltzmann constant, p 
is the mean molecular weight and c s is the isothermal 
sound-speed. Now the evolution is characterized by an 
isothermal shock followed by a weaker, 'D-type' ioniza- 
tion front. The front veloci ty is now i>d < ai on . For a full 
analysis see e.g. IShul (|1991f) . As the hot gas expands, its 
density is reduced. At the same time, the cold surround- 
ing gas is compressed. Under the assumption that the 
homogeneously ionized region consumes all UV-photons 
of the source it follows from Eq. [2] that the density of the 
hot gas for a constant flux and a constant temperature 
Thot at any given time is 



Phot(i) 



m 2 p F Ly 
a B x{t) 



(4) 



To calcul ate the front position x(t) w e follow the ap- 
proach of iDopita fc Sut herland (|2003f ). Under the as- 
sumption of a thin shock which is traveling at a speed v s 
the ram pressure in the hot, ionized gas has to be equal 
to the ram pressure in the cold gas 



where 



Pion -Pold ) 



-Pcold = PO V s = po 



(5) 
(6) 



The pressure on the ionized side of the shock is mainly 
given by the thermal pressure of the hot gas 



-Pion — /-Phot — fPhotC s hot) 



(7) 



where we have already introduced a constant fitting fac- 
tor / to account for the approximations made (e.g. the 
one leading to Eq. |4}. Combining Eq. [6] and [7] and using 
Eqs. [2] and U yields 



x 4 



. dx 
dt 



(8) 



With the initial condition x(to) — x s we can solve by 
integration and obtain 



V 4 x s 
Using Eq. |4]it follows that 

V 4 x s 



(9) 



(10) 



for a plane-parallel infall of a constant flux onto a homo- 
geneous medium. 

2.2. Numerical Method and First Tests 

In order to investigate the effect of different initial 
conditions and levels of UV-radiation on the formation 
of pillars we conduct a parameter study. All simula- 
tions were performed with iVINE (G09a), an implemen- 
tation of ionizing radiation in the tree-SP H code VINE 
(jWetzstein et al.ll2009t iNelson et al.ll2009ft . As the mean 
free path of the atoms and electrons in the cold and the 
hot phase is of ord er A ~ 10 12 cm and A ~ 10 14 cm, re- 
spectively (see e.g. IShul (|1991l) ). which is much smaller 
than the lenght-scales involved in our simulations, the 
gas can be treated by the means of fluid dynamics. The 
evolution of a turbulent molecular cloud under the influ- 
ence of ionization spans several orders of magnitude in 
density. We therefore chose to solve the hydrodynamic 
equations with the method SPH, a Lagrangian approach 
with adaptive resolution. T o prevent artificial fragmenta- 
tion (|Bate &: BurkertHl997f ) the Jeans length has always 
to be resolved with at least 50 particles. This leads to a 
resolution limit which we ensure to be small enough by 
using a sufficient amount of total particles (see M3.7[l . 

In iVINE, the ionizing radiation is assumed to impinge 
plane-parallel onto the simulated volume from the nega- 
tive x-direction. From the surface of infall the radiation 
is propagated along the x-direction by a ray-shooting al- 
gorithm. Along these rays the ionization degree 771 is 
calculated for each particle i. According to the ioniza- 
tion degree the pressure Pj of the particle is calculated 
by a linear interpolation between the temperature Thot 
of the hot, ionized and the temperature T co id of the cold, 
un-ionized gas. Here, we assume both gas components to 
be isothermal, since for the density range considered in 
our simulations heating and cooling should b alance each 
other to approximate isothermality (see e.g. IScalo et al.1 
I1998D . Following Eq. [3] the new pressure in our simula- 
tion is given as 



Pi 



Ti on r]i T nion (l - ?7i)\ k B Pi 



Pu 



Pu 



(11) 



where p\ is the SPH-density of the particle i and p lon — 
0.5 and p n i on = 1-0 are the mean molecular weights of 
the ionized and the un-ionized gas in the case of pure 
hydrogen, respectively. 

As a first test we verify Eq. |9] and fit a value for / by 
ionizing a slab of atomic hydrogen with a constant homo- 
geneous density of p co id = 300 m p cm -3 and temperature 
of Tcoid = 10 K. We perform three different runs, corre- 
sponding to a low flux (-FLy = 1-66 x 10 9 7cm -2 s _1 ), an 
intermediate flux (F^y — 5 x 10 9 7cm _2 s _1 ) and a high 
flux (^y = 1-5 x 10 10 7 cm" 2 s" 1 ). This corresponds 
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t [kyr] 

Fig. 1. — Front position versus time for the three test simulations 
with a different flux impinging on a homogeneous medium. Green, 
blue and red line: simulations with a low, intermediate and high 
flux, respectively. Black lines: solution according to Eq. [9] dotted 

/ = 1, dashed / = 

to the ionization penetrating immediately into the first 
0.55%, 1.67% and 5% of the region, respectively (see Eq. 
[2]). At the same time this is equal to placing the simula- 
tion volume further away or closer to the source, e.g. the 
O-star. The simulations are conducted with the same ac- 
curacy and setup as in the parameter study given below 
(see ij2.3p . Fig. [T] shows the resulting evolution of the 
front. As one can clearly see the approximations lead- 
ing to Eq. [2] (/ = 1, the dotted lines in Fig. Q} do not 
produce satisfactory results. Instead, assuming 

pu = yf-FW (12) 

(i.e. / = the dashed lines in FigJTJ) perfectly 

matches the simulations during the entire simulated time 
of £ s im — 500 kyr. Thus, we keep this assumption for this 
work. 

2.3. Initial Conditions 

To produce different turbulent initial conditions we use 
the same approach as in G09b. We set up 2 x 10 6 par- 
ticles to resemble a homogeneous medium with /o co id = 
300 m p cm' 3 in a volume of (4pc) 3 . Then, we add a su- 
personic velocity field (Mach 10 in most cases) with a 
steep power-law E(k) oc fc~ 2 on the modes k = 1..4. Be- 
fore switching on the ionizing source each setup decays 
freely under the influence of isothermal hydrodynamics 
and periodic boundary conditions until the desired ini- 
tial Mach number is reached (after t s=s 0.8 — 1.0 Myr, de- 
pending on the specific simulation). This self-consistent 
evolution of the turbulence leads to a combination of 
solenoidal and compressible modes, which seems to be 
the case in turbulent molecular clouds 1 . The individ- 
ual particle time-steps in iVINE are determined as in 
G09b by using an accuracy parameter of r acc = 1.0 and 
a Courant-Friedrichs-Lewy (CFL) tolerance parameter of 
tcfl = 0.3. An additional time-step criterion based on 
the maximum allowed change of the smoothing length 
with an accuracy parameter of Th = 0.15 is also em- 
ployed. 

For a more detailed i nvestigation of the role of compressional 
and solenoidal modes see IFederrath eTaLl ll2010h 



After the desired Mach number is reached, the ioniza- 
tion is turned on. Now, the boundaries are still periodic 
in the y- and z-direction. In the negative x-direction the 
boundary is reflecting to represent conservation of flux 
towards the star, whereas in the positive x-direction the 
gas is allowed to stream away freely. To test this ap- 
proach we perform one simulation with open boundaries 
in both x-directions. Gravitational forces are calculated 
without boundaries. This is valid as the free-fall time of 
the whole simulated area is tg s=s 3 Myr, which is much 
longer than the simulation time of tfinai = 0.5 Myr. For 
the tree-based calculation of gravitationa l forces we use 
a mu lti-pole acceptance criterion (MAC, ISpringel et al.l 
1 2 If ) with a tree accuracy parameter of 9 = 5 x 10 -4 . 
The correct treatment of the ionization and the resulting 
acceleration of the particles is guaranteed by the modified 
CFL-condition discussed in G09a. The recombination 
of the hot gas is modeled assuming «b = 2.59 x 10~ 13 
and the cross-section for the ionizing photons is set to 
a = 3.52 x 10~ 18 cm 2 . The initial properties of the dif- 
ferent simulations are listed in Table [TJ The simulations 
were performed on a SGI Altix 3700 Bx2 supercomputer, 
the calculation of each setup took approximately 100 wall 
clock hours on 16 CPUs. 

3. RESULTS OF THE PARAMETER STUDY 
3.1. General Properties 

The time evolution of the density for turbulent regions 
with Mach numbers of 1.5, 5, 7 and 12.5, respectively, is 
shown in Fig. [5] In all of the simulations forming struc- 
tures the same effect as described in G09b takes place. 
The ionizing radiation penetrates deeper into the tur- 
bulent cloud along low-density channels. As the ionized 
gas reacts to its increase in pressure it starts to compress 
the adjacent, un-ionized, higher density regions, thereby 
widening the channels of low density and thus allowing 
the ionization to penetrate even further. We call this 
process 'radiative round-up'. At t = 250 kyr the pre- 
existing high-density structures have been enhanced by 
the outside compression and pillars start to become vis- 
ible. After t = 500 kyr the pillars have achieved charac- 
teristic shapes which match the observations remarkably 
well. Fig. [2] (row 3) and Fig. [3] show the density pro- 
jected along the z-axis at this stage for all simulations of 
the parameter study. 

For a quantitative discussion we investigate the most 
prominent pillar structure in each simulation in more 
detail. To define the tip we take the particle closest 
to the source of radiation 2 above a threshold density of 
Ptresh = 10 4 77ipCm~ 3 . We then take its surrounding, the 
region spanning 1 pc in the x-direction and .3 pc in the 
negative and positive y- and z- direction. The cold, un- 
ionized gas in this region is defined as the pillar 3 . This 
definition allows us to extract the important quantities 
of the most prominent structure by the same algorithm 
for all simulations. The characteristic values, which al- 
low for a comparison with the observations as well as a 
deduction of the underlying physics, are given in Tabled 
for the defined pillar at tfi na i = 500 kyr. Denote that due 

2 In low density, Mach 7 and small box the second tip was taken 
to produce comparable results. 

3 Since we only take the unionized gas inside the surrounding the 
actual volume of the pillar changes from simulation to simulation. 
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a In the case Mach 12.5 we start with an initial velocity fiel d of 
Mach 20 which is then allowed to decay freely as prescribed in i|2.3l 

'To allow for comparable results the box is 2pc in the y- and z- 
direction but 4pc in the x-direction. Therefore, the particle number 
is increased 

TABLE 1 

Listing of the different initial conditions. Given are initial mass, average density and size of the simulation. In addition, 
the impinging flux, turbulent mach number, the largest driving mode of the turbulence and the temperature are listed. 

Mach 5 (G09b) is the standard case as presented in G09b. 
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4.5 ± 1.5 


0.10 


41.8 


2.42 


Mach 4 


12.6 


3.62 


1.32 


20.4 


1.3 ±0.6 


3.9 ±0.9 


0.13 


39.8 


2.20 


Mach 5 (G09b ) 250 kyr 


17.0 


3.60 


1.43 


20.5 


1.8 ± 1.7 


4.0 ± 1.7 


0.16 


60.6 


3.46 


Mach 5 (G09b) 


12.6 


4.56 


1.52 


20.5 


1.1 ±0.7 


4.8 ±0.9 


0.12 


43.4 


1.90 


Mach 7 


14.0 


5.36 


1.71 


20.6 


1.1 ±0.5 


4.0 ±0.9 


0.11 


39.2 


1.46 


Mach 12.5 


8.1 


5.66 


2.56 


20.7 


1.3 ± 1.9 


4.0 ± 1.7 


0.09 


43.5 


1.54 


Mach 5 warm 


10.2 


0.37 


0.23 


19.7 


2.0 ±0.8 


3.3 ± 1.7 




43.6 


2.36 


M5 warm high flux 250 kyr 


7.1 


2.59 


1.0 


20.3 


3.5 ±3.6 


7.0 ±3.7 


0.11 


116 


0.90 


M5 warm high flux 


20.5 


2.38 


0.84 


20.2 


2.3 ± 1.7 


8.0 ±2.1 


0.22 


92.0 


0.77 


low flux 


12.4 


3.75 


1.38 


20.5 


0.9 ±0.5 


3.3 ±0.7 


0.13 


26.5 


1.41 


high flux 250 kyr 


10.9 


6.73 


2.80 


20.8 


2.4 ±3.3 


6.3 ±3.3 


0.09 


108 


3.22 


high flux 


10.9 


7.23 


2.72 


20.8 


2.3 ± 1.3 


8.2 ±2.2 


0.09 


73.5 


2.02 


low density 250 kyr 


3.1 


4.46 


1.82 


20.6 


3.2 ±3.8 


7.0 ±4.0 


0.06 


38.7 


1.74 


low density 


5.79 


2.80 


0.99 


20.3 


2.0 ±0.9 


6.1 ± 1.9 


0.10 


19.0 


1.36 


smaller fcm ax 


2.78 


3.44 


1.16 


20.4 


1.2 ±0.7 


6.2 ±0.9 


0.06 


42.2 


2.45 


small box 


3.39 


2.76 


1.00 


20.3 


1.6 ±0.8 


3.4 ± 1.5 


0.08 


38.2 


2.77 


big box 


51.1 


4.76 


1.94 


20.6 


1.0 ±0.5 


2.9 ±0.7 


0.24 


42.5 


1.78 



TABLE 2 

Results of the parameter study at t = 500 kyr. Listed are the mass, mean density, mean surface density, corresponding 

COLUMN DEPTH, VELOCITY DISPERSION AND THE X- VELOCITY AWAY FROM THE SOURCE OF THE MOST PROMINENT STRUCTURE. THEN THE 
MEAN DIAMETER (SEE EQ. 1131 ) OF THE PILLAR AND THE MEAN DENSITY OF THE HOT GAS ARE GIVEN. FINALLY, THE PRESSURE DIFFERENCE 
AP (SEE EQ. [14)| IS LISTED. FOR THE FIDUCIAL SIMULATION AS WELL AS THE MORE RAPIDLY EVOLVING SIMULATIONS THESE QUANTITIES 

ARE GIVEN AT AN EARLIER STAGE (t = 250 kyr) AS WELL. 



to the adaptive nature of SPH we have a very high res- 
olution in the prominent structures. The pillar in Mach 
5 (G09b) e.g. contains 5.6 x 10 4 particles, which corre- 
sponds roughly to a spatial resolution of 178 x 17.8 x 17.8 
on the 1 pc x 0.1 pc x 0.1 pc of this pillar. This resolu- 
tion is up to now unprecedented and allows for a detailed 
comparison of the kinematics of this pillars with the ob- 
servations (see ^5|). 



We calculate the diameter of a pillar via 

/ M 

^pillar = 2\ (13) 

W pXTT 

with x = 1 pc as the length of the pillar. In addition, 
the mean density 4 of the hot gas is listed in Table [5] It 
is notable that all simulations with the same impinging 

4 Denote that we always give the real density p of the hot gas, 
not the number density n to avoid the factor of p = 0.5 when 
comparing the low density, un-ionized gas to the ionized gas. 
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-5.6 



-4.J 



-4.2 



-3.5 



-2.8 



Mach 1 .5 
: t = kyr 





IMP 




Fig. 2. — Time evolution of four different initial conditions. We show the density projected along the z-axis of four different simulations 
at subsequent stages. From left to right column by increasing Mach number, as indicated. Colour-coded is the surface density, each figure 
is 4pc X 4pc. The significant structures only form above a certain level of turbulence (M > 2) and get less stable with increasing Mach 
number. 



flux share a similar density of the pillar as well as of the 
hot gas. 

As both gas phases are treated isothermal (cf Eq. [3J 
the pressure difference at tfi na \ = 500 kyr is: 

a Pion, final 2Xi on Pion, final f-i a\ 

Affinal = — ■ = — . (14) 

± pillar, final nion Ppillar, final 

APflnai is very close to unity for all simulations 5 . There- 
fore, we derive as a first result that the pillars are in 
thermal equilibrium with the hot surrounding gas. 

3.2. Resolution and Boundary Conditions 

The first two simulations do not address different phys- 
ical properties but rather numerical details. Simula- 
tion low resolution was performed with exactly the same 
setup as Mach 5 ( G09b ), but with eight times less par- 
ticles. This leads to a two times lower spatial resolu- 
tion. Nevertheless, the morphology is comparable to the 

5 In fact, the value is always slightly above one, but this can be 
e.g. attributed to the com plete neglection of the turbulent motion 
in the cold gas in Eq. 1141 



high resolution case (Fig [3] panel 2). The only notice- 
able difference is that the second largest structure in this 
case has already merged with the third structure. Fur- 
thermore, tiny structures are less frequent. The phys- 
ical properties (see Table [2j are similar as well. The 
structures in the low resolution case tend to be a bit 
more massive, which can be expected. Altogether, the 
morphology and the global physical properties are com- 
parable and thus, we conclude that Mach 5 (G09b) is 
reasonably converged. 

In the other test case we investigate the boundary con- 
dition in the negative x-direction. In open boundaries 
this boundary is not reflecting. Instead the gas is al- 
lowed to stream away freely. This leads subsequently to 
a lower density in the ionized region. As a consequence 
open boundaries (Fig. G2 panel 3) looks similar to high 
flux (panel 8, see EJ3T5J) , as the radiation is able to pene- 
trate further into the computational domain. Neverthe- 
less, the formation of pillars still takes place and is not 
strongly affected. Even density and mass assembly of 
the most prominent structure are alike. Therefore, the 



Mach 5 warm 
t = 500 kyr 





Fig. 3. — The final stage of the different simulations. Color-coded is the surface density of the simulations given in Table[2]as indicated 
in the top left corner, each figure is 4pc X 4pc. Denote that in panel 6, 8 and 9 an earlier stage is depicted since the evolution of the 
simulation is more rapid.. Since all panels have the same physical size panel 12 shows only the top left quarter of that simulation which 
includes the most significant structure. 



choice of the boundary condition does not influence the 
overall scenario significantly. As it is more realistic to 
assume hot gas is already present in the region between 
the ionizing source and the simulated part of the molec- 
ular cloud we keep the reflecting boundary condition in 
all other simulations. These reflection can be interpreted 
as flux conservation at the left border of the simulation: 
as much gas streams from the area towards the source 
into the region as is streaming outwards. 

3.3. Turbulent Mach Number 

A main purpose of this study is to disentangle the ef- 
fects of the initial turbulent density distribution on the 
formation of the pillars. Therefore, the level of turbu- 
lence is changed. We take different turbulent setups: 
One has been evolved from a very high Mach number 
(Mach 20) to Mach 12.5. The other four represent dif- 
ferent stages of the decay starting from our fiducial tur- 
bulent setup (G09b) at Mach 10. They are taken at Mach 
7, Mach 5 (G09b), Mach 4 and Mach 1.5, respectively. 
When non-driven turbulence decays, most power is lost 



on the large scale modes. This can be seen in Fig. [2] In 
Mach 5 (G09b) (column 2) and Mach 7 (column 3) the 
surface density is clearly dominated by the large modes, 
which form the prominent fingers. In contrast, in Mach 
1.5 (column 1) no significant pillars evolve, since the ini- 
tial density distribution is already too smooth and the 
dominant mode has decayed to far. This trend can al- 
ready be seen in Mach 4 (Fig. EJ panel 4), where the 
structures are less distinct. Mach 12.5 (column 4) is a 
much more violent case. Since there is a lot of power on 
the largest density scale, structures are evolving. How- 
ever, these are already being torn apart at the same time, 
as discussed in !j4] 

Overall, the evolution is mainly dominated by the pres- 
sure differences between hot and cold gas. Compared to 
the increase in the pressure due to the ionization (three 
orders of magnitude) the differences from varying the 
Mach number are small. However, a small trend is visi- 
ble in the average density of the assembled structure (see 
Table [2]). The higher the Mach number, the higher the 
density of the formed structure. That can be directly 
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related to the density of the initial turbulent filament, 
which is jdenser at a higher Mach number. This effect 
is also visible from the first row of Fig. [5J where the ab- 
solute value in the area of highest density is increasing 
from left to right. 

3.4. Temperature and Pressure 

The most striking difference can be seen between Mach 
5 (G09b) (Fig. [3j panel f) and Mach 5 warm (panel 5). 
Both initial conditions are self-similar, both were set up 
with the same initial random seed for the turbulence and 
they are relaxed until their velocities resemble Mach 5 at 
their respective temperature. Consequently, at the time 
the ionization is switched on, their density distribution 
is identical. Since the impinging flux is the same, the 
radiation ionizes the same regions in both cases. Never- 
theless, in Mach 5 warm no evolution of any filamentary 
structure is visible. This leads to the conclusion that the 
pressure balance between the hot, ionized and the cold, 
un-ionized gas plays a crucial role in the formation of 
structures. 

Taking Eq. [3] and the straightforward assumption that 
only regions with a pressure lower than the pressure of 
the hot gas can be compressed gives 

-fnion, initial ^ Pion, initial (!<-*) 

and thus 

2T- 

Pnion < PionTp^- (16) 
^ nion 

If we assume pi on = 100TO p cm -3 in the beginning, as 
the ionization mainly penetrates the lower density re- 
gions, this equation yields p n ion,ioK < 3.6 x 10 5 m p cm~ 3 
and /9nkm, iook < 3.6 x 10 4 TO p cm -3 . The maximum den- 
sity p m ax = 8.8 x 10 4 m p cm" 3 in both simulations lies in 
between these thresholds. Thus, in Mach 5 (G09b) the 
pressure of the ionized gas is high enough to compress 
even the densest structures, whereas in Mach 5 warm 
several regions are able to resist the compression. There- 
fore, the ionized 'valleys' are not expanding significantly 
in the tangential direction, the density is not lowered as 
much and the ionization is not able to reach much fur- 
ther. At the same time the un-ionized 'hills' are less 
compressed, but since they are closer to the front they 
are accelerated more strongly in the x-direction and the 
initial differences in the front position are leveled out. 

In general, the formation of pillars critically depends 
on whether the density contrast between the dense re- 
gions (phigh), which can not be ionized, and the less dense 
regions (pi ow )) which can be ionized, is lower than the 
temperature ratio between ionized and un-ionized gas. 
By defining the density contrast in the initial conditions 
as Api n it = Phigh/piow and taking into account Eq. [15] 
the critical criterion for the formation of pillars can be 
written as: 

Apinit < (17) 
J- nion 

Since stars will only form in compressed regions, e.g. pil- 
lars, this gives an estimate if a region will undergo trig- 
gered star formation. 

To test this condition we used the same warm initial 
conditions but increased the ionizing flux in M5 warm 
high flux. The flux is now at the level of high flux (see 




o r 1 

100 200 300 400 500 

t [kyr] 



Fig. 4. — The change of the mean density in the hot gas dur- 
ing the simulated time. Solid lines: green low flux, blue Mach 
5 (G09b), red high flux an d ye llow low density. Dotted lines: the 
comparison simulations of i|2.2l at the respective flux. Dashed black 

lines: the analytical solution according to Eq. [4] with / = 

.5[) . This is equivalent to increasing the lowest density 
which can still be ionized, i.e. to increasing p\ ow . As a 
result, Api n it is decreased and thus structure formation 
is possible again in the warm case according to Eq. 1 1 71 
In fact, it can be seen in Fig. [3] (panel 6) that indeed 
pillar formation is triggered again. 

Unfortunately there is no straightforward way to define 
the density contrast in a turbulent medium, especially 
since the impinging flux plays a major role in defining 
'high' and 'low' density as seen in M5 warm high flux. 
However, Eq. [T7]already shows, that increasing the mean 
density p while keeping the temperature constant will not 
help to hinder the formation of pillars. This is supported 
by the fact that pillars form in Mach 5 ( G09b ) as well 
as in low density (see $33}. 

3.5. Initial Flux and Density 

In the next test we vary the impinging photon flux. As 
in the simulations performed in H2.21 the flux is able to 
ionize immediately 0.55%, 1.67% and 5% of a medium at 
a constant density of p = 300 m p cm~ 3 in low flux, Mach 
5 ( G09b ) (intermediate flux) and high flux, respectively. 
The evolution of the density in the hot component is 
shown in Fig. |4] Although the medium is highly turbu- 
lent Eq. QJJstill gives a very good estimate of this density, 
especially after an initial phase. Only the case with the 
high flux differs from the analytical solution. This can 
be understood as Eq. |4] depends on both the penetration 
depth and the density of the ionized gas. A higher flux 
can ionize a larger fraction of the computational volume 
straight away. Equivalently, the evolution of high flux 
would follow Eq. 4 more closely if the mean gas density 
would be higher, thus resulting in a shorter penetration 
length. As long as the penetration length is relatively 
small (< 5%), the turbulent case is still comparable with 
the case of a constant flux and Eq. [TU] still represents 
a valid description of the evolution of a turbulent HII 
region. 

From a morphological point of view, the pillars in the 
simulations with a higher flux (i.e. when the computa- 
tional volume is located closer to the ionizing source) are 
smaller than in the case of a lower flux (Fig. [3j panel 8, 
panel 1 and panel 7 in decreasing flux order). In addi- 
tion, they gain more momentum away from the source 
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and move faster away from the source as the photo- 
evaporation is stronger. At the same time the density 
of the hot gas is higher, leading to denser, more com- 
pressed structures with a smaller diameter. Due to the 
higher photo-evaporation rate their average masses are 
as well lower (see Table [2|) . 

Changing the initial flux is expected to have a simi- 
lar effect as changing the mean density, as p oc x 2 (cf 
Eq. [2]). In low density we reduced the mean density 
by a factor of three. At the same time we reduced the 
flux by a factor of three to avoid an extremely high level 
of ionization degree. In total, this corresponds to the 
same penetration length as in high flux. Thus, we ex- 
pect a similar morphology to evolve, but the densities 
should be lower. In Fig. [3] (panel 9) this can be clearly 
seen. The morphology is similar to high flux, the front 
is at a similar position. Again, the density (Fig. rjj in 
the hot gas evolves similarly to the expectation for a ho- 
mogeneous medium. The mass assembled in the most 
prominent structure (Table [2]) is lower and the density 
of the structure fits the findings of pressure equilibrium 
(see g33). 

Combining these findings with the results of §3. II allows 
us to make an interesting prediction. As the density of 
the hot gas behaves similarly to the case of a homoge- 
neous medium and as the structures are in approximate 
pressure equilibrium with the surrounding hot gas, we 
can predict the density of the structures from the ini- 
tial mean density of the medium, the flux of the source, 
and the time since the ignition of the source or the posi- 
tion of the ionization front. The density of the forming 
structures is thus given as 



log 10 S [g/cm" 



lo Sio s [g/™ ! ] 



271 

' £ <- i ion 
Ppillar ~ Ttt Pion 
nion 



2T ir 



-Po 1 



5 c. 



s.hot 



(t - t Q ) 



(18) 

where we used Eq. [TO] x s depends on the initial density 
and the impinging flux (see Eq. [5]). As we expect the 
assumption of pressure equilibrium to hold in the case of 
a point source, the three-dimensional expression taking 
into account geometrical dilution is given by 



2Tj on fit' c s,hot f , , % 
Ppillar ~ 7f7—P0 I 1 + J {t ~ t ) 



(19) 



where i?s is the Strbmgren radius (for a detailed de riva- 
tion of the three-dimensional front position see e.g. iShul 

[mil). 

3.6. Turbulent Scale 

To study the effect of the largest scales of the initial 
turbulence (the turbulent input scales) we compare Mach 
5 ( G09b) with smaller fc max , a run in which we populate 
modes k = 4... 8, instead oi k = 1..4 as usual. The re- 
sulting surface density in the first 2 pc facing the star 
is shown in Fig. [5] Already in the initial conditions 
(left column) a clear difference can be seen. Whereas the 
power on the larger k modes leads to large, distinct struc- 
tures (top panel), power on the smaller modes show al- 
ready a much more diversified density distribution (lower 
panel). After t nna i = 500 kyr (right column) the ioniza- 
tion leads to an enhancement of the pre-existing struc- 
ture. The densest filaments survive, while the other ma- 
terial is swept away by the ionization. In Mach 5 ( G09b ) 
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Fig. 5. — Projected surface density along the x-axis. The pro- 
jected slice is always 2pc thick. Left: t = Okyr and the slice 
starts at the surface facing the O-star, right: 500 kyr and the slice 
is adjusted to encompass the substructures. Top row: only modes 
k = 1 — 4 are populated initially, bottom row: only modes k = 4 — 8 
are populated initially. 



(top panel) this leads to an excavation of the few, but 
bigger structures and thus to the creation of few, but dis- 
tinct pillars. On the other hand, in smaller fc max (bottom 
panel) more structures, but of smaller scales survive, 
which leads to several, but more diffuse structures. 

Together with £ 13.31 this shows, that only a strong 
enough turbulent driving on a large enough driving scale 
leads to the formation of coherent structures as seen in 
observations. As has been shown in H3.2I this is not an 
effect of the resolution. The turbulence is well enough 
resolved to allow for small enough modes to produce 
fuzzy structure in Mach 5 (G09b), but the evolution un- 
der the influence of UV-radiation is dominated by the 
larger modes. 

Another possibility to change the input scale of the tur- 
bulence is to simply increase or decrease the size of the 
simulation domain. In big box the box size is doubled 
to 8pc. Since the particle number is kept constant, this 
leads to a factor of two lower spatial resolution. So the 
resolution in the part of the domain shown in Fig[3] (panel 
12) is comparable to the low resolution case low resolu- 
tion (panel 2). In small box the box size is halved to 2 pc. 
This corresponds to doubling the spatial resolution. The 
domain has a smaller extent in the x-direction as well, 
which we compensate by taking two times the evolved 
turbulent box in the x-direction. This is valid, since the 
initial conditions were evolved with periodic boundary 
conditions. The particle number is thus 4 x 10 6 , twice as 
high as in most other cases, small box is the situation in 
between Mach 5 (G09b) and smaller A: max , as the largest 
mode is 2 pc, which corresponds to a k = 2 mode in a 
4pc domain. 

As the density distribution and therefore the density 
contrast is self-similar in all three cases we expect that 
the same regions of the initial conditions will form the 
dominant structures. Only the size of the encompassed 
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big box 


305.3 


0.86 


4.37 ± 1.4 


-2.83 


-1.84 


1.26 


big box 


353.3 


0.90 


4.16 ±0.7 


-2.60 


-1.86 


1.34 


big box 


403.7 


0.78 


4.53 ± 1.1 


-2.58 


0.29 


0.41 


big box 


469.3 


0.83 


4.99 ±0.7 


-2.06 


-1.88 


1.48 


high flux 


429.5 


0.57 


11.27 ± 1.0 


1.86 


-1.05 


0.64 


Mach 5 (G09b) 493.3 


0.72 


3.92 ±0.8 


0.19 


0.48 


-0.79 


TABLE 3 



Listing of the proto-stellar cores forming in the 
different simulations. given are mass, formation time, 
average speed away from the source and their formation 
position. 



region and therefore the size and mass of the pillars 
should change, if the process is indeed scale free. In 
Fig. [3] all three simulations (panel 1, panel 11, panel 12) 
show a clear sign of the largest k = 1 mode. The size 
of the structures formed is linearly dependent on the ini- 
tial box-length or size of the largest k-mode. The values 
of Table [2] confirm the importance of the largest mode. 
In the assembled structures the estimated diameters are 
roughly a factor of two different and the masses vary by 
a factor of four. This is as expected since the regions 
initially encompassed by the radiation should differ by a 
factor of two in the y- and the z-direction. 

Taking these results on the turbulent scale into ac- 
count, we conclude that the mass and size of the pillars 
is directly dependent on the input scale of the turbu- 
lence, e.g the size of the driving process or the size of 
the pre-existing molecular cloud. On average, the most 
prominent structures in our simulations with an interme- 
diate flux are d p n « j^turb, where Xturb i s the largest 
turbulent input mode 6 . 

3.7. Star Formation 

In several simulations triggered dense regions form 
cores and are driven into gravitational collapse. Since 
star formation is not the main goal of this study we 
do not replace them by sink particles. Instead we re- 
move the particles forming a core from the simulation to 
avoid a considerable slowdown of the calculation. Fol- 
lowing G09a we define a core as all gas with a density 
above p cr it = 10 7 m p cm -3 in the region around the den- 
sity peak and remove the particles representing this core. 
The core formation is not a numerical effect, since the 
resolution limit as given by iBate fc Burkertl (|1997t ) is 
Pnum = 3 x 10 s m p cm -3 in the lowest resolution case. 
We give the simulation, mass 7 , formation time, average 
speed away from the source and positions in Table [3l If 
we assume the cores to be decoupled from the rest of the 
cold gas then their position at the end of the simulation 
can be estimated by these position and velocities. All 
of them are still close to the prominent structures, some 
are traveling further inside the structures, some are lag- 
ging behind and would by now be slightly outside of the 
pillar, closer to the source. 

6 In fact the value for low density from Table l2l does not match 
precisely, but from Fig. |3]the factor of as 4 between Mach 5 ( G09b) 
and low density can be seen. 

7 The masses do not differ significantly, as we do not follow the 
further accretion process and at the moment of formation the cores 
are still similar. 



Fig. 6. — Projected surface density of the pillars in the simulation 
with Mach 7. (a) x-y projection, (b) x-z projection. Colour- 
coded is the surface density, each figure is 1 pc X 1 pc. The density 
enhancements inside the pillars at the left hand side correspond 
directly to caps, which are not shadowed by the leading tip on the 
right hand side and vice versa. 
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TABLE 4 

Properties of the globules and tips in the four 
simulations. the structures are chosen from the final 

STAGE AT t = 500 kyr , THE NUMBERS ARE GIVEN FROM TOP TO 

bottom. Listed are the mass, the initial tangential 

VELOCITY (vq) OF THE PARTICLES THAT ARE GOING TO FORM THE 
STRUCTURE LATER. AS WELL AS THEIR FINAL TANGENTIAL 
VELOCITY (t£ 0() ). 

As it can be expected, in big box, where the compressed 
structures are most massive, star formation is most fre- 
quent and happens earliest. There is an age spread 
present as well - the earlier a core forms the closer to the 
source it will be. The other simulations where cores form 
during the first 500 kyr after the ignition of the O star 
are high flux and Mach 5 ( G09b ). The higher flux leads 
to a higher compression and thus an earlier formation of 
a core in high flux, due to the higher photo-evaporation 
rate this core is also less massive and moves faster com- 
pared to the core in Mach 5 ( G09b ). 

Altogether, triggered star formation is very likely in 
this scenario. The cores form at the center of the struc- 
tures, but since their velocities differ from the velocity 
of their parental structure they can be decoupled and ei- 
ther wander further into the trunk or lag behind. This 
depends on the specific environment and does show no 
correlation with e.g. the initial flux. At the time they 
become observable there might be no clear correlation to 
their birthplace any more. 

4. PILLARS, CAPS AND GLOBULES 

Another interesting feature can be found when look- 
ing at different projections of the simulation Mach 7 
(Fig. [B]) • By comparing both surface densities it becomes 
clear that the density enhancements inside the pillars are 
caps, which are exposed directly to the UV-radiation in 
the other projection. Thus, the small steps and wiggles 
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seen in the observations can be explained directly by the 
sweeping up of smaller caps, which are not shadowed by 
the leading tips of the pillars. 

As a last point it is interesting to examine the stability 
of the forming structures. In Table|4]we show the correla- 
tion between the formation of globules or pillars and the 
initial turbulent velocity perpendicular to the plane of ir- 
radiation of the particles forming these structures. Here 
we determine the mass of a globule by integrating over 
its material with a density higher than n c = 10 4 cm -3 . A 
tip is defined similarly by taking all the material above 
n c in a sphere of R c = 0.025 pc around the local den- 
sity maximum. Table |4] lists the most important features 
seen in the final stage of Fig. [5J In the physical prop- 
erties of the tips and globules a striking correlation can 
be seen. The tips have lower initial and final tangential 
velocities, the globules much higher ones. Thus, a pillar 
can only be assembled when the leading blob is moving 
slow enough perpendicular to the irradiation that a sub- 
structure can survive in its shadow. Therefore, an envi- 
ronment with lower Mach number favors the formation 
of stable, massive structures, like in the simulation with 
Mach 5 ( G09b). Mach 7 represents the intermediate case 
where tips and globules form. The case of Mach 12.5 
corresponds to a rather violent scenario. Here, globule 1 
has already decoupled from the rest of the gas. In addi- 
tion, even the most prominent pillar (Tip 1) has a high 
velocity and represents a transient stage, which soon de- 
couples from its stem and gets disrupted during the next 
100 kyr. In contrast, Tip 2 is an exceptional case - due 
to its very low tangential velocity it allows for a stable 
pillar even in this violent environment 8 . 

Furthermore, it is worth noting that the tangential ve- 
locity is not affected strongly by the UV-radiation. This 
is reasonable, since the surviving structures experience 
tangential compression from all sides, so the net effect 
of the ionization balances out. In general, our simula- 
tions are able to explain the formation of the observed 
low mass globules or globulettes. Especially the similar- 
ity between the structures forming in Mac h 7 and the 
struc tures observed in the Rosette Nebula (|Gahm et al.l 
120071 their Fig. 3) is striking. 



5. COMPARISON TO OBSERVATIONS 

5.1. General Properties 

At firs t, we compare t he de nsity of the hot gas to obser- 
vations. iLefloch et al.l (|2002f) estimate the electron den- 
sity of the HH-region in the Trifid Nebula 9 from OIII as 
n e — 50 cm" 3 . In a fully ionized region this corresponds 
directly to the density given in Table [H since n c — njj = 
p/mp. Thus, the observed value is very similar to the 
density found in all simulations with an intermediate flux 
at t = 500 kyr. The average density of the pillars is 
around 10 4 m p cm~ 3 , depending on the individual simu- 
lation. This is in very good agreement w ith recent results 
from Herschel bv lSchneider et aT] (|2010f ). where they find 
a typical average density of 1.1 x 10 4 m p cm -3 in the cold 

8 Both statements have been verified by continuing this simula- 
tion until t = 600 kyr. 

9 The e xciting source of the Trifid Nebula is HD 164492A an 
07V star (Lynds ct al. 1985), so the UV-flux is comparable to our 
simulations. 



structures in t he Ro sette Nebula 10 . In the Trifid Nebula, 
ILefloch et all (|2002h estimate N(CS) w 1.8 x 10 13 cm" 2 
for the column density of the dense core, which corre- 
sponds to log 10 [7V(H 2 )/ cm' 2 ] « 22.5 with their con- 
version factor. To compute the H2 column density we 
apply a conversion factor of \ — A(H 2 )/£ = 0.35, us- 
ing a hydrogen abundance of A = 0.7 and assuming 
that all hydrogen is molecular at these densities. In 
all cases where structures form the column densities are 
around log 10 [A(H 2 ) cm" 2 ] « 20.5 (Tabled]). As this 
is the averaged surface density of the entire structure 
in our simulations it is two orders of magnitude lower 
than the observed values for the dense cores. In the 
tips of the pillars (see e.g. Fig. [7]) the peak surface 
density is log 10 [£ max /(gcm" 2 )] = —1.2 which leads to 
a column density log 10 [A(H 2 )/ cm" 2 ] = 22.12, which is 
in goo d agreement with observations. I Thompson et all 
(|2002f ) estimate log 10 [A(H 2 )/ cm" 2 ] w 21.3 for the most 
prominent of the pilla rs in M16. In RCW 120 (e.g. 
iDeharveng et a l. 2009) condensation 4 seems to be a 
good candidate for triggered star formation. The peak 
surface dens ity is log ln [A(H 2 )/ cm " 2 ! = 22.15 - 22.49. 
In addition, lUrquhart et al.l (|2009l ) investigate a sample 
of 60 bright rimmed clouds and find the column den- 
sities in cases of triggered star formation (that is in 
the cases with photodissociation regions (PDRs)) to be 
log 10 [A(H 2 )/cm- 2 ] = 20.9 - 22.8. 

As all these observations match our simulations, we 
conclude that the evolution of the density of the hot gas is 
in good approximation given by the estimate in the case 
of a homogeneous medium (Eq. I10p . Furthermore, since 
the densities of the compressed structures are reproduced 
as well, we can assume that the structures are indeed in 
pressure equilibrium with the hot, surrounding gas. This 
provides the opportunity to determine the density of the 
hot gas and the compressed structures directly from the 
initial mean density, the flux from the source and the 
time since the ignition of the source or the position of 
the ionization front, respectively. 

5.2. Velocity Field of a Singular Pillar 

Another property to compare between simulations and 
observations are the details of the velocity distribution in 
the pillars. Since turbulence is a highly complex process 
it is hard to precisely predict the outcome of simula- 
tions. Therefore, we do not set up a simulation to match 
some specific observation but instead start with realistic 
initial conditions and then look for an observed counter- 
part of the outcome of our simulation. In the following 
we analyze the velocity structure of a typical pillar in the 
Mach 5 simulation which has a similar morphology and 
mass as the only trunk with well observed line-of-sight 
velocities, the Dan cing Queen (DQ) trunk in NGC 7822 
(|Gahm et al.l 120061 see Fig. l|. The authors give the 
diameter as <i bs ~ 0.12 — 0.15 pc, the total estimated 
mass from 12 CO is M obs 9.2 M Q . This is similar the 
simulated pillar (d s i m « 0.12 pc, M s j m 12.6 Mq, see 
Table \2§- If we subdivide the pillar into a head and 
a body the mass splits up into M|, ca( j, s ~ 7.2 Mq and 
-Mbody,s ~ 5.3 Mq (compared to Mh ea d,o ~ 5.7 M Q and 

10 The structures are in the 'Extended Ridge' which is roughly 
10 pc away from the several O-stars (04 to 09) of NGC 2244. 
Thus, the situation there is as well comparable to our simulations. 
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Fig. 7. — The most prominent pillar of the simulation with Mach 
5 in its reference frame. Colour-coded is the surface density, each 
figure is 0.4 pc X 0.8 pc. (a) x-y projection, (b) x-z projection. 
The superimposed grid denotes the bins along which the line-of- 
sight (LOS) velocities i n Fig. [9] are taken. In order to match the 
observational beam-size Gahm ct al. (2006) each bin is 0.04 pc X 
0.04 pc. 




Fig. 8. — For a better comparis on, the Danc i ng Q ueen (DQ) 
trunk in NGC 7822 as observed by I Gahm et all (12009 1 is shown. 
Depicted is a 0.8 pc X 0.8 pc subset of their Fig. 1, roughly to scale 
with our Fig. [7] 

-^body.o ~ 3.5 Mq in the DQ trunk). 

To enable a more detailed comparison we impose a grid 
across the head in the y-x plane (see Fig. [7]) to resemble 
the beams along which the line-of-sight (LOS) velocity 
was taken. We divide the LOS-velocity, ranging from 
v z = — 4kms _1 to v z — 4kms _1 into 80 equally sized 
bins. In each of the velocity bins the mass is integrated. 
Fig. [9] shows the profiles obtained in that way. As we 
do not take any radiative transfer, temperature depen- 
dencies or chemistry into account in our profiles, they 
are not as smooth and symmetric as the observed HCO + 
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Fig. 9. — Line-of-sight(LOS) velocities along the most prominent 
pillar in the simulation with Mach 5. The LOS-velocities are taken 
in bins of the size 0.04 pc X 0.04 pc (see Fig. 01 to match the 
observational resolution. T he m ean velocity in each bin is depicted 
by the dotted lines. In Fig. HOI this velocities are plotted for column 
2, 3 and 4. Denote that in addition to the overall patt ern, the line- 
width is in very good agreement with the observations Gahm et al. 
(2006). 
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Fig. 10. — The line-of-sight(LOS) velocity along the pillar. We 
show the LOS velocity of the most prominent pillar in the simula- 
tion with Mach 5 as function of position along three longitudinal 
cuts. The x-axis is parallel to the major axis. Stars, diamonds and 
triangles correspond to cuts parallel to the x-direction through the 
center (diamond) and the left (star) and right (triangle) side. For 
more details, see Figs. [7] and [9] The location of the head is at 
x = Opc. Top panel (a): projection along the y-axis. Bottom (b): 
projection along the z-axis. A velocity gradient along the pillar 
is clearly visible. The gradien t matches the obs ervations of the 
Dancing Queen Trunk very well[Gahm ct al. (2006:). The fact that 
the velocities for the different cuts rarely cross, especially in the 
bottom plot, is a signature of rigid rotation. This, in combination 
with the overall gradient produces the so called cork-screw pattern, 
as observed. 
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Fig. 11. — The main structure in the radiation driven implosion 
(RDI) model (G09a). Colour-coded is the surface density, each 
figure is 0.4 pc X 0.8 pc. (a) x-y projection, (b) x-z projection. The 
superi mpo sed grid denotes the bins along which the LOS- velocities 
in Fig. 1121 are taken. 
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Fig. 12. — LOS-velocities along the main structure in the RDI- 
model (G09a). The LOS-velocities are taken in bins of the size 
0.04 pcx 0.04 pc (see Fig. 1 11 I t to match the observational resolution. 
Neither the overall p attern, nor the line -width are in agreement 
with the observations Galim ct al. (2006). 



profiles. Nevertheless, the similarities are striking. 

In our simulation the standard mean deviation for all 
lines with a peak higher than 5 x 10 -3 M is on av- 
erage 0.38 km s A which corresponds to a FWHM of 
0.94 km s -1 . This is very similar to the observed FWHM 
of about 1.2 km s . Thus, the irregular motions in the 
simulations and observations correspond very well. Our 
simulations are missing the offset of vq rs — 16.4 km s A 
the speed at which the observed DQ is moving with re- 
spect to us. In addition, the profile is reminiscent of a ro- 
tational pattern, as the peak is shifting from left to right, 
as well as a gradient along the x-direction, since the peak 
is shifting from top to bottom. This so called 'corkscrew' 
pattern (see Fig. [TU)) has often been attributed to mag- 
netic fields. However, our pure hydrodynamical simu- 
lations reproduce the pattern, indicating that magnetic 
fields might play a minor role. Instead, the pattern is 
produced by confining the motion of the turbulent cold 
gas inside the pillars by the hot gas surrounding it. 

As elongated, pillar-like structures can also be the re- 
sult of t he ionization of pre-e xisting clumps (RDI-model, 
see e.g. iMackev fc Lim|[2010t ). we applied the same anal- 
ysis to our previous simulation of this scenario (G09a). 
We took the case of a low ionizing flux impinging onto 
a marginally stable Bonnor-Ebert-Sphere. The low flux 
case was taken to allow for more moderate velocities. 
This scenario results in an elongated feature depicted in 
Fig. Q~T] at t = 600 kyr. The corresponding LOS-profilcs 
are given in Fig. [121 Here, the profiles differ significantly 
from the observations as well as from Fig. [9] First of all, 
there is a double peak. This shows that the structure is 
produced by the collision of two shock-fronts, one mov- 
ing away from us and one moving towards us. These two 
shock fronts, which are encompassing the original den- 
sity enhancements, can be seen directly colliding in the 
perpendicular projection (Fig. ITTa). 

Even if the two peaks become indistinguishable, the 
lines are much broader. Furthermore, there is no de- 
tectable density gradient or rotational pattern visible. 
In addit ion, the veil seen b etween the two smaller pillars 
in Ml 6 ()Hester et al.l Il996h po ses a real challenge. E ven 
the sophisticated simulations of lMackev fc LinJ ( 20101) do 
not reproduce it. In our simulations (see e.g. Fig. |6|) veils 
arise naturally due to the turbulent motions. This is a 
strong indication that the pillars in M16 or in NGC 7822 
are produced by the interplay of pre-existing turbulent 
structures and ionizaing radiation. 

6. CONCLUSIONS AND DISCUSSION 

With iVINE, a tree-SPH code including ionization, we 
perform a parameter study on the formation of pillar- 
like structures and triggered star formation in HII re- 
gions. First we show that our simulations are converged 
and the choice of boundary conditions does not affect 
the outcome significantly (£ I3.2|) . After that, we show 
that ionizing radiation imposed on a turbulent molecular 
cloud can result in the formation of pillar like structures 
which resemble observed pillars in size, mass, density as 
well was velocity structure. Especially the rotational pat- 
tern observed in several pillars indicates that these where 
formed by the ionization of a turbulent cloud and not by 
the RDI of preexisting clumps. 

We thus conclude: 



1. Varying the turbulent Mach number between 1.5 — 
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12.5 changes the morphology. If the Mach number 
is too low, there are no pre-existing structures that 
can be enhanced, if it is too high, the structures 
disperse quickly. The most favorable regime for 
the formation of pillars is Mach 4 — 10. However, 
the physical quantities such as the mass assembled 
and the column density are only weekly dependent 
on the Mach number (see £13.31) . 

2. The formation of pillar-like structures in the case 
with an favorable Mach number critically depends 
on the initial density contrast. Structures and 
therefore stars only form if the density contrast is 
lower than the temperature contrast between the 
hot and the cold gas: Ap init < ^P™- (see q3.4[) . 

3. The evolution of the ionized mass, density and the 
front position in a turbulent medium under the in- 
fluence of different initial fluxes is remarkably sim- 
ilar to the evolution in a homogeneous case. Thus, 
the size and density inside a Hll-region solely de- 
pends on the initial flux, density and the time since 
the ignition of a star or the distance from the star 
and globally follows the simplified analytical pre- 
scription (see 33.5[) . 

4. The density of the resulting pillars is determined 
by a pressure equilibrium between the hot and the 
cold gas. Therefore, the expected density of the 
structures can be calculated as well (see q3.1[) . 

5. The size of the evolving structures critically de- 
pends on the driving modes of the turbulence. 
Smaller driving modes lead to smaller structures. 
In our simulations the relation is roughly d p n w 
^driving/40 (see 33U). 

6. Core and star formation is likely to occur. The 
higher the mass in the structures and the higher 
the initial flux, the earlier cores form (see M3 . T[) . 

Combining 3) and 4) allows us directly to determine the 
density of the forming structures as function of the initial 
mean density of the medium, the flux of the source, and 
the time since the ignition of the source or the position 
of the ionization front (see Eq. [T8]) . 

One has to keep in mind that our approach is sim- 
plified. First of all, no scattering is taken into account. 
Once a electron recombines, the emitted photon is as- 
sumed to be absorbed in the direct neighborhood (on 
the spot approximation). Thus, the reheating of shad- 
owed regions by the adjacent hot gas is not taken into 



account. How much this affects the formation of pillars is 
the topic of ongoing research (Ercolano & Gritschneder, 
in prep) . In addition, we focus on atomic hydrogen only, 
which makes it impossible to follow the precise temper- 
ature evolution as well as the photodissociation regions 
(PDRs). On the other hand, our simulations indicate 
that the pillars are in pressure equilibrium with the hot 
gas. Therefore, the PDRs might be transition regions 
comparable to a thin shock layer which is not resolved. 
Furthermore, we do not take magnetic fields into account. 
These might have impli cations on the global shapes of 
the HII region (see e.g. iKrumholz et alJl2007| ). Never- 
theless, we are able to reproduce the cork-screw mor- 
phologies in the pillars which were sometimes attributed 
to magnetic fields (see £|5.2[) . Another aspect we neglect 
are stellar winds. Although there is clear observational 
evidence the winds of a massive star interact wi th the 
surrounding ISM (e.g. IWestmoauette et "all [2010). it is 
up to now still unclear how effectively they affect its sur- 
roundings. From our simulations we would estimate that 
stellar winds are of minor importance, maybe mainly en- 
hancing the shock front as soon as a lower density in the 
hot gas allows for the effective driving of a stellar wind. 

Altogether, our simulations are able to reproduce even 
the detailed fine-structure of the pillars that have been 
observed with high resolution. In addition, we find that 
the observed line of sight profiles allow for a clear dis- 
tinction between the radiation driven implosion scenario 
and the 'radiative round-up' presented here. Current ob- 
servations are in favor of our 'radiative round-up' mech- 
anism. Besides, our simulations give a deeper insight 
on the tight correlation between the parental molecular 
clouds size, density and turbulence and the structures ex- 
cavated by the ionizing radiation. The ionization acts as 
a magnifying glass, revealing the condition of the molec- 
ular cloud previous to the ignition of the massive star. 
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